\part{基本的动态系统}
\chapter{军备竞赛模型}
\section{基本的模型设置}
模型总是包含三个部分：内生变量，外生变量和参数。下面是一个军备竞赛模型。
\begin{equation}\label{eq_ds_arm}
\begin{bmatrix}
	\Delta x_{1t}\\ 	\Delta x_{2t}
\end{bmatrix}=\underbrace{\begin{bmatrix}
		-\alpha& \beta \\ \gamma & -\delta
\end{bmatrix}}_{\bm{A}}\begin{bmatrix}
	x_{1t}\\ 	 x_{2t}
\end{bmatrix}+\underbrace{\begin{bmatrix}
		\theta&0\\ 0 & \eta
\end{bmatrix}}_{\bm{B}}\begin{bmatrix}
	z_{1,t}\\ z_{2t}
\end{bmatrix} 	
\end{equation}


其中，$ x_{1t},x_{2t} $是内生变量，表示两个国家的武器存量，可以看到武器存量的变化依赖于自身现有的武器存量和敌国的武器存量。$ z_{1t},z_{2t} $是两个外生变量。$ \alpha,\beta,\gamma,\delta,\theta,\eta $是参数。

根据附录\ref{chp_de}中差分方程理论，这个系统的稳态为，
\begin{equation}\label{eq_ds_st}
\begin{bmatrix}
	\bar{x}_{1,t}\\ 	\bar{x}_{2,t}
\end{bmatrix}=-\bm{A}^{-1}\bm{Bz}	
\end{equation}


该方程的稳定性，可以通过对矩阵$ \bm{A} $的特征值的讨论而得到。

一个动态模型建立后，通常会进行两个方面的分析：\textbf{一是}脉冲响应分析，即外生冲击对内生变量的影响，通常是一条时间轨迹。我们把它称为脉冲响应分析。\textbf{二是}参数对内生变量的影响，我们把它成为敏感性分析。接下来，我们通过对参数的赋值来看到这一点。
\section{数值解}
我们首先给各参数赋值如下，
\[ \alpha = 0.5, \beta = 0.25,\gamma  =0.25,\delta = 0.5, \theta=1,\eta=1 \]

然后，基于这些赋值，可以得到矩阵$ \bm{A} $的特征值为，
\[ \lambda_{1}=-0.25,\lambda_{2}=-0.75 \]
很明显，$ \lambda_{i}+1<1, i \in 1,2 $，本系统在这些参数的设置下是稳定的。

再次，我们给外生变量赋一个初值，
\[ z_{1,t}=1,z_{2,t}=1 \]

然后，根据\eqref{eq_ds_st}式就可以得到稳态值为，
\[ \bar{x}_{1t}=\bar{x}_{2t}=4 \]

\subsection{脉冲响应分析}
现在让我们改变外生变量，从而产生一个冲击，然后观察该冲击下各内生变量的变化轨迹，即如何从一个稳态抵达另一个稳态。我们改变外生变量$ z_{2t}=2 $，其他不变，可以称此时的外生变量取值为终值。在这样的取值下，系统的特征值不受影响，但稳态发生改变，对应的稳态终值为，
\[ \bar{x}_{1t}=6.67,\bar{x}_{2t}=5.33 \]

我们现在关心的是系统如何从前一个稳态，怎样变化来到这个新的稳态。我们可以利用\eqref{eq_ds_arm}式来模拟内生变量的改变路径。首先把差分项展开，
\begin{equation}\label{eq_ds_arm2}
\bm{x}_t=(\bm{A}+\bm{I})\bm{x}_{t-1}+\bm{Bz}_t	
\end{equation}

这样就有一个比较清晰的动态方程。然后令前一个稳态，即$ \bar{\bm{x}}=(4,4)' $为初值，$ \bm{z}_t=(2,1)' $，根据\eqref{eq_ds_arm2}式逐步往前迭代，就可以得到一系列的$ \bm{x}_t $，基本上20期以后，系统就已经稳定下来了。
\begin{figure}[H]
\includegraphics[scale=0.7]{arm_irf.pdf}
\caption{外生冲击下内生变量的脉冲响应图}
\end{figure}


\subsection{敏感性分析}
\paragraph{全局稳定}现在我们来改变参数，再次观察系统的动态。改变$ \alpha = 0.7 $，此时系统特征值为，
\[ \lambda_{1}=-0.33,\lambda_{2}=-0.87 \]
足见系统仍然是稳定的。

系统的稳态值为，
\[ \bar{x}_{1}=2.61,\bar{x}_{2}=3.3 \]

类似的，我们可以继续考察此时内生变量的变化轨迹。\eqref{eq_ds_arm2}式的初值仍然为前一个稳态值$ (4,4)' $，外生变量不变为$ (1,1)' $，但矩阵$ \bm{A} $已经发生了改变，根据\eqref{eq_ds_arm2}式继续模拟，可以得到下图，
\begin{figure}[H]
	\includegraphics[scale=0.7]{arm_sen1.pdf}
	\caption{参数改变下的内生变量的变化轨迹图}
\end{figure}

\paragraph{鞍点稳定}前面参数的变化，系统仍然处于全局稳定中，现在让参数的变化使得系统处于鞍点稳定中。现在参数按如下取值，
\[ \alpha=0.25,\beta=0.5,\gamma=0.5,\delta=0.25,\theta=\eta=1 \]

外生变量取值为，
\[ z_1=z_2=-1 \]

此时系统的特征值为，
\[ \lambda_{1}=0.25,\lambda_{2}=-0.75 \]
可见其中一个$ \lambda_{1}+1>1,\lambda_{2}+1<1 $，因此是鞍点稳定。系统的稳态经过计算，发现依然是，
\[ \bar{x}_{1}=\bar{x}_{2}=4 \]

因为鞍点稳定的特点是只有内生变量恰好处在鞍点路径上时，才会收敛到稳态，其他位置，内生变量都是发散的。


现在来考虑在一个外生冲击下，内生变量如何脉冲响应。画这个脉冲响应图时，有一个小技巧，鞍点稳定下，当冲击发生后，计算第二期数据时要让其中一个变量迅速发生改变，让其调整到鞍点路径上，这样，两个内生变量才能重新又回到稳态。那么其中那个内生变量要如何改变呢？

要得到鞍点路径，根据附录中的\eqref{de_eq_x}式，可知内生变量的路径可以写为，
\[ \begin{bmatrix}
	x_{1t}\\ x_{2t}
\end{bmatrix} = \bm{v}_1(\lambda_{1}+1)^tc_1+\bm{v}_2(\lambda_{2}+1)^tc_2+\begin{bmatrix}
\bar{x}_1\\\bar{x}_2
\end{bmatrix} =\bm{v}_1(0.25)^tc_1+\bm{v}_2(1.25)^tc_2+\begin{bmatrix}
\bar{x}_1\\\bar{x}_2
\end{bmatrix}\]
其中，$ \bm{v}_1,\bm{v}_2 $是对应$ \lambda_{1},\lambda_{2} $的特征向量，$ c_1,c_2 $是和稳态相联系的两个常数。要得到鞍点路径，必须要去掉一个大于1的特征值。因此，上式对应的鞍点路径为，
\[ \begin{bmatrix}
	x_{1t}\\ x_{2t}
\end{bmatrix} = \bm{v}_1(0.25)^tc_1+\begin{bmatrix}
	\bar{x}_1\\\bar{x}_2
\end{bmatrix}\]
这个方程是不够的，我们需要的是一个把$ t $期和$ t-1 $期联系起来的方程。但上述方程是很有启发的，我们可以把它写成一个一般形式，对于$ t $期和$ t+1 $期，有，
\begin{align}\label{eq_ds_xt}
	 \begin{bmatrix}
		x_{1t}\\ x_{2t}
	\end{bmatrix} =& \bm{v}_1(\lambda_{1}+1)^tc_1+\begin{bmatrix}
		\bar{x}_1\\\bar{x}_2
	\end{bmatrix}\\\nonumber
	 \begin{bmatrix}
	x_{1,t+1}\\ x_{2,t+1}
\end{bmatrix}  =&\bm{v}_1(\lambda_{1}+1)^{t+1}c_1+\begin{bmatrix}
	\bar{x}_1\\\bar{x}_2
\end{bmatrix}	
\end{align}

两式相减意味着，
\begin{equation}\label{eq_ds_x2}
 \begin{bmatrix}
	\Delta x_{1t}\\\Delta x_{2t}
\end{bmatrix} = \lambda_{1}\bm{v}_1(\lambda_{1}+1)^tc_1 	
\end{equation}
另一方面\eqref{eq_ds_xt}式意味着，
\begin{equation}\label{eq_ds_x3}
\begin{bmatrix}
	x_{1t}-\bar{x}_{1t}\\ x_{2t}-\bar{x}_2
\end{bmatrix} = \bm{v}_1(\lambda_{1}+1)^tc_1 	
\end{equation}

联立\eqref{eq_ds_x2}与\eqref{eq_ds_x3}式，有，
\begin{equation}\label{eq_ds_saddle}
 \begin{bmatrix}
	\Delta x_{1t}\\\Delta x_{2t}
\end{bmatrix} = \lambda_{1}\begin{bmatrix}
x_{1t}-\bar{x}_1\\ x_{2t}-\bar{x}_2
\end{bmatrix}	
\end{equation}
这就是内生变量在鞍点路径上的运行轨迹。现在假设外生变量发生一个改变，比如$ z_1 $从-1变到-0.5，内生变量如何在第二期调整到鞍点路径上？首先根据系统的动态方程，有，
\begin{equation}\label{eq_ds_saddle2}
\Delta x_{1t}=-\alpha x_{1,t}+\beta x_{2,t}+\theta z_{1,t}	
\end{equation}

如果仅考虑上式，然后以上一个稳态$ (4,4)' $替代上式的$ (x_{1,t},x_{2,t})' $，那么$ t+1 $期的$ x_{1,t+1} $可能并不一定在鞍点路径上。所谓鞍点路径，本质上就是同期变量$ x_{1,t} $和$ x_{2,t} $间的一个函数关系。该函数关系可以通过把\eqref{eq_ds_saddle}式第一个方程的$\Delta x_{1t} $代入到\eqref{eq_ds_saddle2}式而得到。经过这种替换，有，
\begin{equation}\label{eq_ds_saddle3}
x_{1t}=\frac{\beta}{\alpha+\lambda_{1}}x_{2,t}+\frac{\theta}{\alpha+\lambda_{1}}z_{1t}+\frac{\lambda_{1}}{\alpha+\lambda_{1}}\bar{x}_1 	
\end{equation}

现在，就比较好办。比如基于任何一个点或者说就是上一个稳态$ (4,4) $可以利用\eqref{eq_ds_saddle2}式得到$ x_{1,t+1} $，而该点对应在鞍点路径上$ x_{2,t+1} $可以通过\eqref{eq_ds_saddle3}式来得到。一旦有了在鞍点路径上的点，就可以按照\eqref{eq_ds_saddle}式继续运动即可。

\begin{figure}[H]\centering
\includegraphics[scale=0.3]{鞍点跳跃示意图.jpg}
\caption{鞍点路径跳跃示意图}
\end{figure}

写代码的时候，先依照\eqref{eq_ds_saddle2}式算出下一期内生变量$ x_{1,t+1} $的值，再依照\eqref{eq_ds_saddle3}式将内生变量$ x_{2,t+1} $调整到鞍点路径，再依照\eqref{eq_ds_saddle}式运动系统即可。

\begin{figure}[H]
	\includegraphics[scale=0.7]{arm_irf2.pdf}
	\caption{鞍点路径下的脉冲响应图}
\end{figure}

本章全部代码见\lstinline|arms_race.R|。
\chapter{IS-LM 的动态模型}
\section{模型的基本设置}
一个动态的IS-LM模型可以书写如下：
\begin{align}\label{eq_IS_lm}
m_t-p_t=&\psi y_t-\theta i_t\\\label{eq_IS_is}
 y_t^d=&\beta_0-\beta_1(i_t-\Delta p_t^e)\\\label{eq_IS_p}
\Delta p_t=&\mu (y_t-y_t^n)\\\label{eq_IS_y}
\Delta y_t=&v(y_t^d-y_t)
\end{align}
这里，$ m $是货币总量的对数，$ p $是价格水平的对数，$ y^d $是总需求的对数，$ y $是总供给(或总产出)的对数,$ y^n $是潜在产出水平的对数，$ i $是名义利率。这个模型的内生变量有四个，即$ y,p,y^d,i $。还包含了外生变量三个，即$ \beta_0,m,y^n $。第一个方程是LM曲线，表达了货币总量不变的情况下，总产出与利率的正比关系。第二个方程是IS曲线，表达了总需求和利率的反比关系。第三个方程是菲利普斯曲线，表明价格对产出的反应，$ \mu $是调整的速度。第四个方程表达的是产出对总供给和总需求的差异的调整方式。

现在我们把上述四个方程的模型调整成一个仅包含2个差分方程的系统，即价格$ p $的差分方程和产出$ y $的差分方程。因此，我们要消掉其他两个内生变量$ y^d $和$ i $。这一点可以通过对\eqref{eq_IS_lm}和\eqref{eq_IS_is}式的联立而得到，
\begin{align}\label{eq_IS_yd}
	y^d_t=&\beta_0-\frac{\beta_1\psi}{\theta}y_t+\frac{\beta_1}{\theta}(m_t-p_t)+\beta_1\Delta p_t^e\\
	=&\beta_0-\frac{\beta_1\psi}{\theta}y_t+\frac{\beta_1}{\theta}(m_t-p_t)+\beta_1\Delta p_t\\\label{eq_IS_i}
	i=&-\frac{1}{\theta}(m_t-p_t-\psi y_t)
\end{align}

注意到第二个等号中我们消去了$ \Delta p_t^e $这个期望项。因为期望，我们一般认为是理性预期，一般是通胀加一个白噪声。但我们此处认为是完美预测，无不确定性，因此，$ \Delta p_t=\Delta p_t^e $。将$ y_t^d $的解带入\eqref{eq_IS_y}式，然后关于价格的差分方程\eqref{eq_IS_p}式本身不需要改动，整个系统可以重新写为，
\begin{align*}
\Delta y_t=&v\left[\beta_0-\left(\frac{\beta_1\psi}{\theta}+1\right)y_t+\frac{\beta_1}{\theta}(m_t-p_t)+\beta_1\Delta p_t\right]\\
=&v\left[\beta_0-\left(\frac{\beta_1\psi}{\theta}+1\right)y_t+\frac{\beta_1}{\theta}(m_t-p_t)+\beta_1\mu (y_t-y_t^n)\right]\\
\Delta p_t=&\mu (y_t-y_t^n)
\end{align*}

再次重新整理成我们熟悉的方式，有，
\begin{equation}\label{eq_IS_py}
\begin{bmatrix}
	\Delta p_t\\
	\Delta y_t
\end{bmatrix}=\underbrace{\begin{bmatrix}
		0&\mu\\ \frac{-v\beta_1}{\theta}&v\left(\beta_1\mu-\frac{\beta_1\psi}{\theta}-1\right)
\end{bmatrix}}_A \begin{bmatrix}
	p_t\\ y_t
\end{bmatrix}+\underbrace{\begin{bmatrix}
		0&0&-\mu\\ v&\frac{v\beta_1}{\theta}&-v\beta_1\mu
\end{bmatrix}}_B\begin{bmatrix}
	\beta_0\\ m_t\\ y_t^n
\end{bmatrix}	
\end{equation}


\section{参数设置和稳态分析}
\begin{table}[H]\centering
	\caption{参数和外生变量设置}
\begin{tabular}{llll}\hline
	符号&值&	符号&值\\\hline
$\psi$&0.05&$ \theta $&0.5\\
$ \beta_1 $&50&$ \mu $&0.01\\
$ v $&0.2&&\\\hline
$ m_0 $&100&$\beta_0$&2100\\
$ y_t^n $&2000&&\\\hline
\end{tabular}
\end{table}


在这些参数设置下，\eqref{eq_IS_py}式可以写成，
\[ \begin{bmatrix}
	\Delta p_t\\
	\Delta y_t
\end{bmatrix}=\underbrace{\begin{bmatrix}
		0&0.01\\ -20&-1.1
\end{bmatrix}}_{\bm{A}} \begin{bmatrix}
	p_t\\ y_t
\end{bmatrix}+\underbrace{\begin{bmatrix}
		0&0&-0.01\\ 0.2&20&-0.1
\end{bmatrix}}_{\bm{B}}\underbrace{\begin{bmatrix}
	2100\\ 100\\ 2000
\end{bmatrix}}_{\bm{z}_t} \]

首先根据附录\ref{chp_de}的稳态公式，可以得到系统的稳态为，
\[ \begin{bmatrix}
	\bar{p}_t\\ \bar{y}_t
\end{bmatrix}=-\bm{A}^{-1}\bm{Bz}_t=\begin{bmatrix}
1\\ 2000
\end{bmatrix} \]

当你有了这些稳态以后，你可以把它们带入\eqref{eq_IS_yd}和\eqref{eq_IS_i}式得到其他内生变量$ i_t,y_t^d $的稳态为，
\[ \bar{i}=2,\qquad \bar{y}^d=2000 \]

然后，我们依然可以计算$ \bm{A} $的特征值，
\[ \lambda_{1}=-0.23,\lambda_{2}=-0.87 \]
很明显$ \lambda_{1}+1<1,\lambda_{2}+1<1 $，我们此时的系统是全局稳定的。

\section{脉冲响应和敏感性分析}
按照前面分析的步骤，可以观察，当货币供给$ m_0 $从100增加到101时，系统如何从当前稳态变动到另一个稳态。
\begin{figure}[H]
\includegraphics[scale=0.7]{m_irf.pdf}
\caption{$ m_0 $的脉冲下各内生变量的响应}
\end{figure}

当产出对货币需求的弹性$ \psi $从0.05下降到0.01时，系统如何当前稳态变动到另一个稳态，
\begin{figure}[H]
	\includegraphics[scale=0.7]{m_sen.pdf}
	\caption{$ \psi $变化下各内生变量的变动轨迹}
\end{figure}

\paragraph{代码实现}\lstinline|is_lm.R|
\begin{lstlisting}[language = R, mathescape=false]
	rm(list = ls())
	library(patchwork)
	library(tidyverse)
	library(ggplot2)
	library(matlab)
	# parameters setting
	psi <- 0.05
	beta1 <- 50
	v <- 0.2
	theta <- 0.5
	mu <- 0.01
	
	# exo variables
	mt <- 100
	beta0 <- 2100
	ytn <- 2000
	
	# 常微分方程求解结果封装
	ode_solve <- function(psi, beta1,v, theta, mu, mt, beta0,ytn){
		# coef matrix
		A <- matrix(c(0,mu,-v*beta1/theta,v*(beta1*mu-beta1*psi/theta-1)),2,byrow = T)
		B <- matrix(c(0,0,-mu,v,v*beta1/theta,-v*beta1*mu),2,byrow =T)
		zt <- matrix(c(beta0,mt,ytn),ncol = 1)
		
		# steady state
		st <- -solve(A) %*% B %*% zt
		# eig value
		ev <- eigen(A)$values
		return(list(A=A, B=B, zt=zt, st=st, ev=ev))
	}
	
	rlt <- ode_solve(psi, beta1,v, theta, mu, mt, beta0,ytn)
	rlt
	
	# IRF: m0
	mt <- 101
	new_rlt <- ode_solve(psi, beta1,v, theta, mu, mt, beta0,ytn)
	
	x <- matrix(0,nrow = 4, ncol = 20, dimnames = list(c('p','y','yd','i')))
	for (i in 1:20) {
		if ( i == 1){
			x[1:2,1] <- rlt$st
		}else {
		# xt = (A+I)x_t-1+ Bz
			x[1:2,i] <- (new_rlt$A + eye(2)) %*% x[1:2,i-1] + new_rlt$B %*% new_rlt$zt 
		}
		x[4,i] <- -(mt-x[1,i]-psi*x[2,i])/theta
		x[3,i] <- beta0-beta1*(x[4,i]-mu*(x[2,i]-ytn))
	}
	picdata <- t(x) %>% as.data.frame()
	p1 <- ggplot(picdata, aes(x = 1:nrow(picdata),y = y)) + geom_line() + labs(x='') +theme_bw()
	p2 <- ggplot(picdata, aes(x = 1:nrow(picdata),y = p)) + geom_line() + labs(x='')+ theme_bw()
	p3 <- ggplot(picdata, aes(x = 1:nrow(picdata),y = yd)) + geom_line() + labs(x='')+ theme_bw()
	p4 <- ggplot(picdata, aes(x = 1:nrow(picdata),y = i)) + geom_line() + labs(x='')+ theme_bw()
	(p1 + p3)/(p2+p4)
	ggsave('../m_irf.pdf')
	
	# sensitivity:
	psi <-  0.01
	sen_rlt <- ode_solve(psi, beta1,v, theta, mu, mt, beta0,ytn)
	
	x <- matrix(0,nrow = 4, ncol = 30, dimnames = list(c('p','y','yd','i')))
	for (i in 1:30) {
		if ( i == 1){
			x[1:2,1] <- rlt$st
		}else {
			x[1:2,i] <- (sen_rlt$A + eye(2)) %*% x[1:2,i-1] + sen_rlt$B %*% new_rlt$zt
		}
		x[4,i] <- -(mt-x[1,i]-psi*x[2,i])/theta
		x[3,i] <- beta0-beta1*(x[4,i]-mu*(x[2,i]-ytn))
	}
	picdata <- t(x) %>% as.data.frame()
	p1 <- ggplot(picdata, aes(x = 1:nrow(picdata),y = y)) + geom_line() + labs(x='') +theme_bw()
	p2 <- ggplot(picdata, aes(x = 1:nrow(picdata),y = p)) + geom_line() + labs(x='')+ theme_bw()
	p3 <- ggplot(picdata, aes(x = 1:nrow(picdata),y = yd)) + geom_line() + labs(x='')+ theme_bw()
	p4 <- ggplot(picdata, aes(x = 1:nrow(picdata),y = i)) + geom_line() + labs(x='')+ theme_bw()
	(p1 + p3)/(p2+p4)
	ggsave('../m_sen.pdf')	
\end{lstlisting}
\chapter{汇率超调}
\section{理论分析和方程设置}
在开放经济体中，有两个国际均衡关系：
\begin{itemize}
	\item 货币市场由非抛补利率平价决定，即本国利率高于外国利率，市场预期本币未来将贬值。名义汇率是完全弹性的，它伸缩自如。名义汇率是直接标价法，汇率上升，本币贬值。
	\item 商品市场由购买力平价决定，即名义汇率等于国内价格和国外价格之差。但价格调整是刚性的，购买力平价只在长期成立。
\end{itemize}

\begin{definition}{汇率超调的描述}{exc}
当货币供给增加时，汇率会立刻到达一个比它新的稳态更高的水平，然后随着价格调整到新的稳态，汇率逐渐下降了。
	\end{definition}

本质上，货币供给增加会立刻使得名义利率下降，在利率平价下，名义汇率会有升值压力。但由购买力平价决定的商品市场是由国内外价格水平决定，这个价格水平即期没有变化，因此，要满足汇率升值预期，当期汇率必须立马贬值，也就是说汇率立马上升了。但由于价格的刚性，这个上升是过度的。随着价格开始往上缓慢调整，新的价格水平下，上期的汇率相对此时就过高了，于是随着价格的调整，汇率逐渐恢复至应有的均衡水平。这个均衡水平比原来最初的均衡水平仍然是要高的（购买力平价）。


经济的动态可如下描述：
\begin{align*}
	m_t-p_t=&\psi y_t-\theta i_t\\
	y_t^d=&\beta_0+\beta_1(s_t-p_t+p_t^*)-\beta_2i_t\\
	\Delta p_t=&\mu(y_t-y_t^n)\\
	\Delta s_t^e=&i_t-i_t^*
\end{align*}
其中，$ m,p,y $分别是货币供给、价格水平和产出水平的对数，$ i $是名义利率，$ y^d $是需求水平的对数，$ s $是汇率的对数，$ y^n $是潜在产出的对数，$ p^*,i^* $是外国的价格水平和名义利率，$ s^e $是汇率的期望。这里，我们假设$ s^e=s $。

类似的，第一个方程是LM曲线，这个和上一章的封闭经济一样。第二个方程是关于产品市场的IS曲线，该方程发生了变化，因为需要考虑外国的需求的变化。该需求正向依赖于真实汇率(由购买力平价的差异所定义)，反映了以价格水平体现的外部经济的竞争性。该需求还负向依赖于名义利率。第三个方程是菲利普斯曲线，第四个方程是利率平价。

这个系统有五个内生变量需要决定：$ p_t,y_t,i_t,y_t^d,s_t $。但我们只有4个方程，我们需要额外增加一个方程来描述经济。从模型来看，主要是$ y_t $的行为如何决定没有详细设定，要么给它一个设定，要么直接排除掉这个内生变量。为了简化，我们通过设置$ y_t=y_t^n $来直接排除该变量。但这样的操作会使得第三个方程中价格的变化是0，这不是我们想要的，因此，在价格调整方程中设置$ y_t^d=y_t^n $看来是最佳的。这样，经济的方程结构可以书写如下，

\begin{align*}
	m_t-p_t=&\psi y^n_t-\theta i_t\\
	y_t^d=&\beta_0+\beta_1(s_t-p_t+p_t^*)-\beta_2i_t\\
	\Delta p_t=&\mu(y^d_t-y_t^n)\\
	\Delta s_t^e=&i_t-i_t^*
\end{align*}

现在，我们需要重新排列上述方程组，来获得价格水平和名义汇率的差分方程。即要在这两个差分方程中消掉另外两个内生变量$ i $和$ y^d $。这可以通过联立第一和第二个方程得到，
\begin{align}\label{eq_exc_i}
	i_t=&-\frac{1}{\theta}(m_t-p_t-\psi y_t^n)\\\label{eq_exc_y}
	y_t^d=&\beta_0+\beta_1(s_t-p_t+p_t^*)+\frac{\beta_2}{\theta}(m_t-p_t-\psi y_t^n)
\end{align}

然后将其代入第三和第四个方程，可以最终得到如下方程，
\begin{equation}\label{eq_exc_ps}
\begin{bmatrix}
	\Delta p_t\\ \Delta s_t
\end{bmatrix}=\underbrace{\begin{bmatrix}
		-\mu(\beta_1+\frac{\beta_2}{\theta})&\mu\beta_1\\
		\frac{1}{\theta}&0
\end{bmatrix}}_{\bm{A}}\begin{bmatrix}
	p_t\\ s_t
\end{bmatrix}+\underbrace{\begin{bmatrix}
		\mu & \frac{\mu\beta_2}{\theta}&-\mu\left(\frac{\psi\beta_2}{\theta}+1\right)&0&\mu\beta_1\\
		0&-\frac{1}{\theta}&\frac{\psi}{\theta}&-1&0
\end{bmatrix}}_{\bm{B}}\begin{bmatrix}
	\beta_0\\ m_t\\ y_t^n\\ i_t^*\\ p_t^*
\end{bmatrix}	
\end{equation}


一旦获得上式，那么接下来的分析套路是显然的。
\section{校准和稳态}
\begin{table}[H]\centering
	\caption{参数校准和外生变量赋值}
\begin{tabular}{lll|lll}\hline
	参数&含义&值&外生变量&含义&值\\\hline
$ \psi $& 货币需求弹性& 0.05&$ m_0 $&货币供给&100\\
$ \theta $&利率的半弹性&0.5&$ \beta_0 $&加总需求的截距&500\\
$ \beta_1 $&$ y_t^d $对真实汇率的弹性&20&$ y_t^n $&潜在产出&2000\\
$ \beta_2 $&$ y_t^d $对利率的弹性&0.1&$p_0^*$&国外价格水平&0\\
$ \mu $&价格调整速度&0.01&$ i_0^* $&国外名义利率&3\\\hline
\end{tabular}
\end{table}

根据稳态公式，
\[ \begin{bmatrix}
	\bar{p}_t\\ \bar{s}_t
\end{bmatrix}=-\bm{A}^{-1}\bm{Bz}_t =\begin{bmatrix}
1.5 \\ 76.51
\end{bmatrix}\]

\paragraph{稳定性分析}$ \bm{A} $的特征值为-0.74，0.53，因此$ \bm{A}+\bm{I} $的特征值为0.26和1.54，这是一个鞍点稳定。因此，对于鞍点稳定而言，一定要存在一个完全弹性的变量或者说跳跃变量，它可以瞬间抵达鞍点路径，从而整个系统收敛。而这个跳跃变量很明显就是汇率。

\section{脉冲响应分析和敏感性分析}
现在假设货币供给从100增加到101，看看系统如何变化。此时，根据\eqref{eq_exc_ps}式，从旧稳态运动到下一期的值为，
\[ p_{t+1}=1.5,s_{t+1}=74.51 \]

这一对值一般不在鞍点路径上，必须要迅速调整汇率使得其在鞍点路径上。根据\eqref{eq_ds_saddle}式，对于鞍点路径必然满足，
\begin{equation}\label{eq_exc_sadd}
\begin{bmatrix}
	\Delta s_t\\ \Delta p_t
\end{bmatrix}=\lambda_{1}\begin{bmatrix}
	s_t-\bar{s}\\ p_t-\bar{p}
\end{bmatrix}	
\end{equation}

其中，$ \lambda_{1} $是绝对值小于1的那个特征值。将上式中关于$ s_t $的第一个方程与\eqref{eq_exc_ps}式第二个方程联立就能得到鞍点路径上$ p_t$和$s_t $的函数关系为，
\[ s_t=-\frac{-\lambda_1\theta\bar{s}+i_t^*\theta-\psi*y_t^n+m_t-p_t}{\lambda_{1}\theta}\]
从而可以计算得到$ s_t=80.20 $。即这一期新的点应该是$ p_t=1.5,s_t=80.20 $，然后再依据\eqref{eq_exc_sadd}式运动，这样，我们就可以得到关于$ s_t,p_t $的运动轨迹。再利用\eqref{eq_exc_i}和\eqref{eq_exc_y}式，就可以得到另外两个内生变量的轨迹，我们将其绘制如下图所示，
\begin{figure}[H]
\includegraphics[scale=0.7]{exc.pdf}
\caption{货币供给冲击下的脉冲响应图}
\end{figure}

从中可以看到汇率首先一个跃迁，然后缓慢回落至一个比原来高一点的水平上。

然后我们改变参数，重新观察脉冲响应。我们把价格的调整速度$ \mu $从0.01调整到0.001，让价格调整更为缓慢，重新进行前面的脉冲响应分析，图如下， 从中可以看到汇率的调整更加缓慢了。
\begin{figure}[H]
	\includegraphics[scale=0.7]{exc_sen.pdf}
	\caption{更换慢的调整速度下货币供给冲击后的脉冲响应图}
\end{figure}
 
 上述图形实现的R代码如下所示(\lstinline|exch_overshoot.R|)。
 \begin{lstlisting}[language = R, mathescape=false]
rm(list = ls())
library(ggplot2)
library(matlab)
library(patchwork)
# parameters setting
psi <- 0.05
theta <- 0.5
beta1 <- 20
beta2 <- 0.1
mu <- 0.01

# exo variables
m0 <- 100
beta0 <- 500
ytn <- 2000
p0star <- 0
i0star <- 3

ode_solve <- function(psi, theta, beta1,beta2, mu, m0, beta0,ytn,p0star,i0star){
	# coef matrix
	A <- matrix(c(-mu*(beta1 + beta2/theta),mu*beta1,1/theta,0),2,2,byrow = T)
	B <- matrix(c(mu,mu*beta2/theta,-mu*(psi*beta2/theta+1),0,mu*beta1,0,-1/theta,psi/theta,-1,0),nrow = 2,byrow = T)
	z <- matrix(c(beta0,m0,ytn,i0star,p0star),ncol = 1)
	
	# steady state
	st <- -solve(A) %*% B %*% z
	# eig value
	ev <- eigen(A)$values
	return(list(A=A, B=B, z=z, st=st, ev=ev))
}

# 稳态
rlt <- ode_solve(psi, theta, beta1,beta2, mu, m0, beta0,ytn,p0star,i0star)
rlt$st

# 特征值
rlt$ev

# 脉冲响应分析
# 新的稳态
m0 <- 101
rlt_new <- ode_solve(psi, theta, beta1,beta2, mu, m0, beta0,ytn,p0star,i0star)
st_new <- rlt_new$st

# 冲击下，第一期的值
jmp <- (rlt_new$A + eye(2)) %*% rlt$st + rlt_new$B %*% rlt_new$z
# 汇率的瞬间跳跃
jmp[2] <- -(-rlt_new$ev[1]*theta*st_new[2]+i0star*theta-psi*ytn+m0-jmp[1])/(rlt_new$ev[1]*theta)

# 鞍点路径的运动
picdata <- data.frame(p = c(rlt$st[1],jmp[1]),s = c(rlt$st[2],jmp[2]))
for (i in 3:15) {
	picdata[i,] <- 1/(1-rlt_new$ev[1])*picdata[i-1,]-rlt_new$ev[1]/(1-rlt_new$ev[1])*st_new
}
picdata$i <- -(m0-picdata$p-psi*ytn)/theta
picdata$i[1] <- -(100-picdata$p[1]-psi*ytn)/theta # 修正初期值,即用原来的m0

picdata$ytd <- beta0 + beta1*(picdata$s-picdata$p+p0star)+beta2/theta*(m0-picdata$p-psi*ytn)
picdata$ytd[1] <- beta0 + beta1*(picdata$s[1]-picdata$p[1]+p0star)+beta2/theta*(100-picdata$p[1]-psi*ytn)
picdata$time <- 1:nrow(picdata)

pexc <- ggplot(picdata, aes(x = time, y = s)) + geom_line() + theme_bw()
pp <- ggplot(picdata, aes(x = time, y = p)) + geom_line() + theme_bw()
py <- ggplot(picdata, aes(x = time, y = ytd)) + geom_line() + theme_bw()
pi <- ggplot(picdata, aes(x = time, y = i)) + geom_line() + theme_bw()
(pexc / pp)|(py/pi)
# ggsave('../exc.pdf')


# 敏感性分析
mu <- 0.001
rlt <- ode_solve(psi, theta, beta1,beta2, mu, m0=100, beta0,ytn,p0star,i0star)

# 冲击下新的响应路径
rlt_new <- ode_solve(psi, theta, beta1,beta2, mu, m0=101, beta0,ytn,p0star,i0star)
st_new <- rlt_new$st

# 冲击下，第一期的值
jmp <- (rlt_new$A + eye(2)) %*% rlt$st + rlt_new$B %*% rlt_new$z
# 汇率的瞬间跳跃
jmp[2] <- -(-rlt_new$ev[1]*theta*st_new[2]+i0star*theta-psi*ytn+m0-jmp[1])/(rlt_new$ev[1]*theta)

# 鞍点路径的运动
picdata <- data.frame(p = c(rlt$st[1],jmp[1]),s = c(rlt$st[2],jmp[2]))
for (i in 3:15) {
	picdata[i,] <- 1/(1-rlt_new$ev[1])*picdata[i-1,]-rlt_new$ev[1]/(1-rlt_new$ev[1])*st_new
}
picdata$i <- -(m0-picdata$p-psi*ytn)/theta
picdata$i[1] <- -(100-picdata$p[1]-psi*ytn)/theta # 修正初期值,即用原来的m0

picdata$ytd <- beta0 + beta1*(picdata$s-picdata$p+p0star)+beta2/theta*(m0-picdata$p-psi*ytn)
picdata$ytd[1] <- beta0 + beta1*(picdata$s[1]-picdata$p[1]+p0star)+beta2/theta*(100-picdata$p[1]-psi*ytn)
picdata$time <- 1:nrow(picdata)

pexc <- ggplot(picdata, aes(x = time, y = s)) + geom_line() + theme_bw()
pp <- ggplot(picdata, aes(x = time, y = p)) + geom_line() + theme_bw()
py <- ggplot(picdata, aes(x = time, y = ytd)) + geom_line() + theme_bw()
pi <- ggplot(picdata, aes(x = time, y = i)) + geom_line() + theme_bw()
(pexc / pp)|(py/pi)
# `ggsave('../exc_sen.pdf')
 \end{lstlisting}

现在，给大家抛出一个有意思的问题：在经济学中，一个典型的家庭将按如下规则决策消费，
\begin{align*}
	\max_{\{C_t\}_{t=0}^T}&\sum_{t=0}^T\beta^t\ln C_t\\
	s.t.\qquad &C_t+B_t=(1+R)B_{t-1}+W_t\\
&	B_{-1}=0\\
& B_T=0
\end{align*}
可以得到一阶条件如下，
\begin{align*}
	C_{t+1}=&\beta(1+R)C_t\\
	B_{t+1}=&(1+R)B_{t}+W_{t+1}-C_{t+1}\\
	=&(1+R)B_{t}+W_{t+1}-\beta(1+R)C_{t}
\end{align*}

写成矩阵：
\[ \begin{bmatrix}
	C_{t+1}\\ B_{t+1}
\end{bmatrix}= \underbrace{\begin{bmatrix}
\beta(1+R)&0\\
-\beta(1+R)&1+R
\end{bmatrix}}_{A+I}\begin{bmatrix}
C_{t}\\ B_{t}
\end{bmatrix}+\begin{bmatrix}
0\\ W_{t+1}
\end{bmatrix}\]

这里会发现如果不满足$ 1= \beta(1+R)$，则系统没有稳态解，除非$ C=0 $，此时才有稳态解。特别地，如果我们取值$ \beta=0.97,R=0.02, W=10 $，会发现这是个鞍点稳定系统，因为矩阵$ A+I $的特征值为1.02和0.98，稳态为$ \bar{C}=0,\bar{B}=-500 $。

此时，我们可能更加感兴趣在某些边界条件下，比如$ B_{-1}=0,B_T=0 $，最优消费的时间路径会是怎样的，而这可能不一定是在稳态附近展开的一条路径。

如果$ T=30 $，上述问题本质上是求解$ \{C_1,C_{2},\cdots,C_{30},C_{31}\} $和$ \{B_1,B_2,\cdots, B_{29}\} $。对于$ t=1 $，一阶条件表明(此处把关于$ B_t $的一阶条件的时间脚标往前推一期)，
\begin{align*}
	C_2=\beta(1+R)C_1\\
	B_1=W_1-C_1	
\end{align*}
对于$ t=2,\cdots,29 $，一阶条件为，
\begin{align*}
	C_{t+1}=&\beta(1+R)C_t\\
	B_{t}=&(1+R)B_{t-1}+W_{t}-C_{t}
\end{align*}
对于$ t=30 $，一阶条件为，
\begin{align*}
	C_{31}=&\beta(1+R)C_{30}\\
0=&(1+R)B_{29}+W_{t}-C_{30}
\end{align*}

因此，一共有方程60个，未知数60个，解此（非）线性方程组即可。
\begin{lstlisting}[language = R, mathescape=false]
rm(list = ls())
library(pacman)
p_load(bvpSolve,ggplot2,magrittr,signal,rootSolve,matlab)
devtools::load_all()

# 参数设置
beta <- .97
R <-  .02
WL <- 10
numT <- 31 # 31期

find_c <- function(cc, parms = list(numT = 30, beta = 0.97, R = 0.02, WL = 10)){
	eps <- B <- numeric(parms$numT)
	# 初值：B0 = 0
	B[1] <- parms$WL - cc[1]
	# B迭代
	for (i in 2:parms$numT) {
		B[i] <- (1+parms$R)*B[i-1] + parms$WL - cc[i]
	}
	
	# 消费迭代
	for (i in 1:(parms$numT - 1)) {
		eps[i] <- cc[i+1] - parms$beta * (1+parms$R) * cc[i]
	}
	# 终值：B(T) = 0
	eps[parms$numT] <- (1+R)*B[parms$numT-1] + parms$WL - cc[parms$numT]
	return(eps)
}

# 解多元非线性方程组
cc <- 10*rep(1,numT)
ans <- multiroot(find_c, start = cc, parms = list(numT = numT, beta = beta, R = R, WL = WL))
picdata <- data.frame(ts = 1:length(ans$root), cc = ans$root)
ggplot(picdata, aes(x = ts, y = cc)) + geom_line()
\end{lstlisting}